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ABSTRACT 

Context. Radiative transfer plays a key role in the star formation process. Due to a high computational cost, radiation-hydrodynamics 
simulations performed up to now have mainly been carried out in the grey approximation. In recent years, multi-frequency radiation- 
hydrodynamics models have started to emerge, in an attempt to better account for the large variations of opacities as a function of 
frequency. 

Aims. We wish to develop an efficient multigroup algorithm for the adaptive mesh refinement code RAMSES which is suited to heavy 
proto-stellar collapse calculations. 

Methods. Due to prohibitive timestep constraints of an explicit radiative transfer method, we constructed a time-implicit solver based 
on a stabilised bi-conjugate gradient algorithm, and implemented it in RAMSES under the flux-limited diffusion approximation. 
Results. We present a series of tests which demonstrate the high performance of our scheme in dealing with frequency-dependent 
radiation-hydrodynamic flows. We also present a preliminary simulation of a three-dimensional proto-stellar collapse using 20 fre¬ 
quency groups. Differences between grey and multigroup results are briefly discussed, and the large amount of information this new 
method brings us is also illustrated. 

Conclusions. We have implemented a multigroup flux-limited diffusion algorithm in the RAMSES code. The method performed well 
against standard radiation-hydrodynamics tests, and was also shown to be ripe for exploitation in the computational star formation 
context. 
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1. Introduction 

Numerical studies of star formation are very demanding, as 
many physical mechanisms need to be taken into account (hy¬ 
drodynamics, gravity, magnetic fields, radiative transfer, chem¬ 
istry, etc). Models are rapidly increasing in complexity, provid¬ 
ing ever-more realistic interpretations of todays highly advanced 
observations of protostellar systems, such as the recent ground¬ 
breaking images taken by the ALMA interferometer^. Radiative 
transfer plays a key role in star formation; it acts as a conduit to 
remove compressional heating during the initial stages of cloud 
collapse, enabling an isothermal contraction (Larson 1969; Ma- 
sunaga et al. 1998), and it also inhibits cloud fragmentation in 
large-scale simulations (see Bate 2012, for instance). State-of- 
the-art simulations thus require the solutions to the full radia¬ 
tion magneto-hydrodynamics (RMHD) system of equations, and 
three-dimensional simulations have only just recently become 
possible with modem computers (see Commer 9 on et al. 201 la; 
Tomida et al. 2013, for example). In particular, including fre¬ 
quency dependent radiative transfer is essential to properly take 
into account the strong variations of the interstellar gas and dust 
opacities as a function of frequency (see for example Ossenkopf 
& Henning 1994; Li & Draine 2001; Draine 2003; Semenov 
et al. 2003; Eerguson et al. 2005). Three-dimensional calcula- 


^ http: //WWW .almaobservatory.org/en/press-room/press-releases/ 
771-revolutionary-alma-image-reveals-planetary-genesis 


tions including full frequency-dependent radiative transfer are 
still out of reach of current computer architectures. 

In order to overcome this difficulty, much effort has been 
spent in recent years developing mathematically less compli¬ 
cated, yet accurate approximations to the equations of radiative 
transfer. Such representations include diffusion approximations, 
the Ml model, short and long characteristics methods, and Vari¬ 
able Eddington Tensor descriptions. One of the simplest and 
most widely used is the flux-limited diffusion (ELD) approxima¬ 
tion (Levermore & Pomraning 1981), and has been applied to 
many areas of physics and astrophysics. 

These methods tend to drastically reduce computational cost, 
but still they are often integrated over all frequencies (also 
known as the ‘grey’ approximation) as the multi-frequency for¬ 
malism remains too expensive. Only in recent years have multi¬ 
group methods (whereby the frequency-dependent quantities are 
binned into a finite number of groups and averaged over the 
group extents) appeared in numerical codes (Shestakov & Offner 
2008; van der Holst et al. 2011; Vaytet et al. 2011; Davis et al. 
2012; Zhang et al. 2013; Vaytet et al. 2013b, to mention a few), 
as a first step towards accounting for the frequency dependence 
of gas and dust quantities in calculations. Multigroup methods 
are very suited to astrophysics as they allow adaptive widths of 
groups, thus enabling to use a small number of groups, con¬ 
centrating frequency resolution where it is needed (i.e. mostly 
where opacity gradients are strong). Group boundaries are usu¬ 
ally chosen once at the start of the simulations, but more complex 
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schemes have also been written with moving adaptive borders 
(Williams 2005), as absorption and emission coefficients gener¬ 
ally vary with the material temperature and density. 

Commer 9 on et al. (2011b) implemented frequency- 
integrated FLD in the adaptive mesh refinement (AMR) code 
RAMSES (Teyssier 2002; Fromang et al. 2006), and devised 
an adaptive time-stepping scheme in a follow-up paper (Com- 
mergon et al. 2014). In this third paper, we extend the method to 
a multigroup formalism. Even though multigroup effects were 
found to only have a small impact on one-dimensional (ID) 
simulations of star formation (Vaytet et al. 2012, 2013a), they 
are expected to be enhanced in 3D where the optical thickness 
can markedly vary along different lines of sight (see Kuiper et al. 
2011). Using a frequency-dependent method is also the only 
way to correctly model ultra-violet radiation from stars being 
absorbed by surrounding dust and re-emitted in the infrared. 
Finally, multigroup formalisms, compared to grey methods, 
are known to significantly affect the structures of radiative 
shocks (Vaytet et al. 2013b), alter energy transport in stellar 
atmospheres (Chiavassa et al. 2011), as well as being essential 
to neutrino transport in core collapse supemovae explosions 
(see Mezzacappa et al. 1998, for instance). 

We first present the numerical method we have used, fol¬ 
lowed by a series of tests against analytical solutions, and we 
end with an application to the collapse of a gravitationally un¬ 
stable cloud, comparing grey and multigroup results. 


2. Numerical method 

The FLD multigroup radiation hydrodynamics equations in the 
frame comoving with the fluid are 


+ V • Ipu] 

dtipu) -h V • Ipu 0 u -r PI] 
dtEr -r V • [u(Pt + P)] 


dtEg + V . [uEg] 


0 

Ng 

Ng 


8=^ 

-rV- 


pKRg 


-VP, 


-Pg : Vu + V 


cA„ 


-VE, 


+KY,gpc{@g{T)-Eg) 

r^g+112 

: I (9y(yPy)Jy 

n 


-rVu 


( 1 ) 


where c is the speed of light, p, u, P, and T are the gas density, 
velocity, pressure, and temperature, respectively. Pt is the total 
energy Ej = pe + Ijlpu^ + Eg (e is the internal specific 
energy), and I is the identity matrix. We also define 
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where X = E,¥ which represent the radiative energy and pres¬ 
sure inside each group g which holds frequencies between yg- 1/2 
and y^+i/ 2 - Efg is the total number of groups, 0^(T) is the en¬ 
ergy of the photons having a Planck distribution at temperature 
T inside a given group. The coefficients Ag, K^g and /cr^ are, re¬ 
spectively, the flux-limiter, the Planck and the Rosseland means 
of the spectral opacity Ky inside a given group. 


We employ a commonly used operator splitting scheme 
whereby the equations of hydrodynamics are first solved explic¬ 
itly using the second order Godunov method of RAMSES includ¬ 
ing the radiative terms involving V • u and VP^, while the evo¬ 
lution of the radiative energy density and its coupling to the gas 
internal energy is solved implicitly (for more details on the dif¬ 
ferent equations which are solved explicitly and implicitly, we 
refer the reader to the exhaustive description in Commergon et al. 
2011b). 

To discretize the equations solved in the implicit step, we 
linearize the source term following 

= ©giT") + - T"), (3) 

where the prime denotes the derivative with respect to tempera¬ 
ture. This then enables us to write a set of discretised equations, 
expressed here in ID for simplicity, for the evolution of the gas 
temperature (where Cy is the gas heat capacity at constant vol¬ 
ume) 
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The terms with superscripts n+\ refer to the variables evaluated 
at the end of the timestep Ar, while superscripts n indicate the 
state at the beginning of the timestep. Subscripts i represent the 
grid cell, and / +1 /2 are for cell interfaces. In addition, V, S , and 
Ax are the cell volume, the interface surface area, and the cell 
width, respectively. The subscripts a and p denote the frequency 
groups where the subscript g is already in use. 

The implicit step requires the inversion of a large matrix 
which holds the system of equations (4) and (5), and this is per¬ 
formed using a parallel iterative method^. In the case of the grey 
approach (and in our specific case of Cartesian coordinates), the 
matrix to invert in the implicit step is symmetric, which allows to 
use the conjugate gradient method (for further details, see Com- 
mer 9 on et al. 201 lb, 2014). However, in the multigroup case, the 
interaction between radiative groups adds non-symmetric terms, 
for which we had to implement a similar but more advanced sta¬ 
bilised bi-conjugate gradient (BiCGSTAB) algorithm (van der 
Vorst 1992). 

Finally, the term which depends on a derivative with respect 
to frequency {dy) accounts for energy exchanges between neigh¬ 
bouring groups due to Doppler effects. It is computed using the 
method described in Vaytet et al. (2011), and treated explicitly 
as part of the final line in equation (11) in Commer 9 on et al. 
(2011b) (note that we now have one such line per frequency 
group). 

^ The use of a direct inversion method is not suited to the very large 
systems of equations that need to be solved in heavy 3D simulations. 
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3. Method validation 

We present in this section the numerical tests performed to assess 
the accuracy of our method. 


3.1. Dirac diffusion 


We consider the one-dimensional two-group radiation diffusion 
equation in a static medium with no coupling to the gas. The 
equations to solve are then 

= » , 6 , 

For a constant pK^ coefficient and a Dirac amplitude value of 
Eq at xq as initial condition, the analytical solution Ea in a p- 
dimensional space is 


Ea{.^-> 0 


Eq 
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where = c/(3p/^R). 

We choose a box of length L = \ cm, where xq = 0.5 cm. 
The gas has a uniform density p = 1 g cm"^. The initial total ra¬ 
diative energy is set to 1 erg cm“^ except in the center (inside the 
two central cells) where the peak value Eq is set to 10^ erg cm“^. 
The two frequency groups’ boundaries are (in Hz) [10^, 10^^] 
and [10^^, 10^^], chosen so that, in the central region, the ra¬ 
diative energy in the first group is about two orders of magni¬ 
tude lower than in the second one. The Rosseland opacity in 
the first group is set to a:ri = 1 cm^ g“^ and 10 times higher in 
the second group. The domain is initially divided into 32 cells 
(coarse grid level of 5) and 4 additional AMR levels are enabled 
(effective resolution of 512). The refinement criterion is based 
on the gradient of the total radiative energy. There is no flux- 
limiter, i.e. Xg = 1/3 for both groups, and a fixed timestep of 

= 2.5 X 10“^^ s is used for the coarse level (since we are using 
the adaptive time-stepping scheme of Commer 9 on et al. 2014, 
the timestep is divided by two per AMR level increase). 

Figure 1 (top) shows the radiative energy density profiles at 
time t = 2 X 10“^^ s. The two radiative energies are in good 
agreement with the analytical curves. The relative error (bottom 
panel) on the total radiative energy is always less than 10%. The 
zones where the relative error on group 2 exceeds the 10% mark 
correspond to regions where this group does not contribute to the 
total energy, making this error largely unimportant. 


3.2. Radiating plane 

Graziani (2008) proposed an analytical multigroup test in spher¬ 
ical geometry. Let’s consider a sphere of radius R, and temper¬ 
ature Ts which is surrounded by a cold medium at temperature 
To < Ts, heat capacity Cy, and spectral absorption coefficient 
pKy = CTy. The test consists in computing the time-dependent 
spectrum at a given distance r > R from the sphere center in the 
absence of any hydrodynamic motion. In the case where the heat 
capacity tends to infinity, the gas temperature is constant, and the 
spectrum is the superposition of the black-body spectrum of the 
cold medium and the spectrum of the radiation emitted by the 
hot sphere. The analytical solution for the spectral energy is 

Ey = B(y, To) + - [B(y, T,) - B(y, To)] E(y, r - Rj) (8) 

r 



Distance (cm) 


Fig. 1. Diffusion test: numerical (circles and squares) and analytical 
solutions (solid lines) at time t - 2x 10"^^ s. Bottom panel: relative 
error. The dotted line corresponds to the AMR levels in the simulation. 


where 



As our grid is Cartesian, we adapted this test to a slab geom¬ 
etry. Instead of a radiating sphere, we consider a radiating plane. 
The analytical solution can then be found making r and R tend to 
infinity, while keeping the distance r-R constant (Gentile 2008). 
We then simply have 

Ey = B(y, To) + [B(y, T,) - B(y, To)] E(y, r-Rj). (10) 

In our simulation, the temperature of the hot slab was set to T^ = 
1500 eV and the medium was at Tq = 50 eV. The domain size 
is 0.1347368 cm with the hot slab located at the left boundary. 
The domain was divided into 32 identical cells. We considered 
60 groups logarithmically evenly spaced in the range [0.5 eV, 
306 keV], a fixed timestep of At = 10“^^ s was used, and no 
flux limiter (i.e. Ag = 1/3) was applied. The gas absorption 

coefficient was set to = 2 x 10^^(p^)“^ cm“^with hyg the 
energy in eV of the middle of each radiative group. Figure 2 
shows the spectral radiative energies obtained compared to the 
analytical ones at a time t = 10“^^ s, sampled at a distance x = 
0.04 cm (this implies that r - R = 0.04 cm in the analytical 
solution), which corresponds to the center of the tenth cell. The 
numerical results are in excellent agreement with the analytical 
solution. 

3.3. Non-equilibrium radiative transfer with picket fence 
model 

Su & Olson (1999) developed analytical solutions for a ID prob¬ 
lem involving non-equilibrium radiative transfer with two radia¬ 
tive energy groups (see also Zhang et al. 2013). Radiative energy 
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Fig. 2. Radiating plane test: numerical (circles) and analytical solu¬ 
tions at time t = 10"^^ s. 


is injected inside a small region of a uniform domain, diffuses 
and heats up the gas. The two groups have different opacities, so 
that the radiative energies propagate at different speeds through 
the medium. There are several assumptions in their analytical 
study. The heat capacity at constant volume is assumed to be 
Cv = deldT = aT^ where a is a parameter. The group inte¬ 
grated Planck distribution is assumed to be 


Bg=Pg 





( 11 ) 


where is the radiation constant, and pg are parameters which 
verify the condition Pg = 1. They then define the dimen¬ 
sionless coordinate x = az where z is the coordinate in physical 
units, and d" = TjgPgCrg. The absorption coefficients CTg are in¬ 
dependent of frequency, and scattering is ignored. The dimen¬ 
sionless time is 


4(2rC(T \ 

Of / ’ 
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and the dimensionless radiative energy density and internal en¬ 
ergy are 

respectively, where To is a reference temperature. 

The radiation source is applied for a finite period of time 
(0 < T < To) inside the region |x| < vq, and gas dynamics are 
neglected. The equations solved are then 


dtEj 


dtEg - V . 
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Following Zhang et al. (2013), we performed the case C of 
Su & Olson (1999) and compared the results with their analytical 



0.1 1 10 100 


Distance x 

Fig. 3. Non-equilibrium radiative transfer test of Su & Olson (1999). 
Profiles from the numerical simulation of Ui (top), U 2 (middle) and 
V (bottom) are shown for t = 3 (solid lines) and t = 30 (dashed 
lines). The results are compared to the analytical solutions of Su & 
Olson (1999) for t = 3 (circles) and t = 30 (diamonds). 


solution for the radiation diffusion. In case C there are two radia¬ 
tion groups, with pi = p 2 = 1/2. The absorption coefficients are 
chosen as cri = 2/101 cm“^ and 0-2 = 200/101 cm“^ and the 
parameter a used to evaluate the heat capacity is a = 4aR. The 
reference temperature is set to Tq = 10^ K. The radiation source 
parameters are vq = 1/2 and tq = 10. To avoid spurious bound¬ 
ary condition effects, we used a computational domain twice the 
size of Zhang et al. (2013), spanning -102.4 < x < 102.4, di¬ 
vided into 2048 identical cells. The left and right boundary con¬ 
ditions were both set to periodic. The initial state of the physical 
variables were p = 1 g cm“^ and T = I K, and matter and radi¬ 
ation were in equilibrium. A fixed timestep At = 0.1 was used, 
and no flux limiter (i.e. Ag = 1/3) was applied. The results 
are shown in Fig. 3, where an excellent agreement between the 
numerical and analytical solutions can be seen. 

3.4. Radiative shocks with non-equilibrium diffusion 

The fourth test looks at the coupling between the fluid mo¬ 
tion and the radiative transfer, solving the complete radiation- 
hydrodynamics equations (Eq. 1). Lowrie & Edwards (2008) 
developed semi-analytic solutions to a ID non-equilibrium dif¬ 
fusion problem involving radiative shocks of various Mach num¬ 
bers M. We simulated the At = 2 and M = 5 cases using 
6 frequency groups. For both runs, the domain was split in two 
uniform regions where the hydrodynamic and radiation variables 
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Fig. 4. Radiative shocks with non-equilibrium diffusion. The left column are the results for the A1 = 2 case, while the right column is the Af = 5 
setup. In all panels the numerical results are marked by symbols (circles and squares), the analytical solutions are represented by solid lines, and 
the mesh AMR level by dotted black lines. Panels (a) and (c) display the gas density as a function of distance. Panels (b) and (d) show the gas 
(black) and total radiative (magenta) temperatures. They also present the radiative temperatures inside the individual groups 6g with more colours 
(see legend in panel d). 


satisfy the Rankine-Hugoniot jump relations for a radiating fluid 
in an optically thick medium. The fluid variables inside the ghost 
cells at the left and right boundary conditions were kept as the 
initial pre- and post-shock state values throughout the simulation 
(imposed boundary condition). The gas has an ideal equation of 
state with a mean atomic weight ju = 1 and a speciflc heat ra¬ 
tio y = 5/3, and the Planck and Rosseland opacities are set to 
pKp = 3.93 X 10“^ cm“^ andpAfR = 0.848902 cm“\ respectively, 
in all frequency groups. We used the HLL Riemann solver for 
the hydrodynamics with a CFL factor of 0.5 and no flux-limiter 
for the radiation solver (i.e. Ag = 1/3). 

For the At = 2 case, the initial conditions in the left 
(pre-shock) region are pl = 5.45887 x 10“^^ g cm"^, t/R = 
2.3547 X 10^ cm s“\ Tr = 100 K and in the right (post-shock) 
region pr = 1.2479 x 10“^^ gcm“^, mr = 1.03 x 10^ cm s“\ 
and Tr = 207.757 K. The domain ranges from -1000 cm to 
1000 cm. Five frequency groups were evenly (linearly) dis¬ 
tributed between 0 and 2 x 10^^ Hz, and the sixth group ranged 
from 2 X 10^^ Hz to inflnity. The domain was initially divided 
in 32 cells and 4 additional AMR levels were enabled (effec¬ 
tive resolution of 512). The refinement criterion was based on 
gas density and total radiative energy gradients. The density and 
temperature (gas and radiation) profiles are displayed in Figs 4a 


and 4b. The thermodynamic quantities (gas density and tem¬ 
perature) are represented by the black symbols, while the radi¬ 
ation temperatures, defined by Og = (Eg/aR)^^^, are marked by 
coloured squares. The numerical solutions show an excellent 
agreement with the semi-analytical solutions of Lowrie & Ed¬ 
wards (2008) (solid lines). The simulation data were shifted by 
13.67 cm to place the density discontinuity at x = 0 for compar¬ 
ison with the analytical solutions; this corresponds to the shift 
the shock suffers as the radiative precursor develops, until the 
stationary state is reached. 

In the At = 5 case, the initial conditions are pl = 5.45887 x 
10“^^ gcm“^, Ur = 5.8868 x 10^ cm s“\ Tr = 100 K, and 
Pr = 1.96405 X 10“^^ gcm“^, ur = 1.63 x 10^ cm s“\ and 
Tr = 855.72 K. The domain ranges from -4000 cm to 4000 cm. 
Five frequency groups were evenly (linearly) distributed be¬ 
tween 0 and 10^"^ Hz, and the sixth group ranged from 10^"^ Hz 
to inflnity. The domain was initially divided in 32 cells and 
7 additional AMR levels were enabled (effective resolution of 
4096). The refinement criterion was based on gas density, gas 
temperature, and total radiative energy gradients. The results in 
Figs 4c and 4d show again an excellent agreement with the semi- 
analytical solutions. The simulation data were this time shifted 
by 193.5 cm to bring the density discontinuity to x = 0. 
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Fig. 5. (a) Strong scaling results for the multigroup diffusion test (cir¬ 
cles) and the ideal MHD Orszag-Tang vortex calculation (diamonds). 
The ideal theoretical speedup is at the boundary between the grey and 
white areas. The Occigen machine has 12-core processors; the 12-core 
region is marked by the green hatched area, (b) Weak scaling efficiency 
for the RHD and MHD runs, (c) CPU time per BiCGSTAB iteration 
for a ID diffusion problem (circles) and a ID non-stationary radiative 
shock (squares) as a function of the number of frequency groups Ng 
used. The curve for a time proportional to is shown with a dotted 
line. 


4. Algorithm performance 

We present in this section some tests to assess the scaling perfor¬ 
mance of our algorithm. 


4.1. Strong and weak scaling 

The strong and weak scaling runs were performed on the 
CINES Occigen^ supercomputer, which uses Intel® E5-2690 
(2.60 GHz) processors. We compared the scaling of our radia¬ 
tive transfer scheme to the native MHD scheme in RAMSES. The 
setups used for the RHD and MHD runs were a 2D version of the 
Dirac diffusion test using 4 frequency groups and a 2D Orszag- 
Tang vortex simulation (Orszag & Tang 1979), respectively. In 
the RHD calculation, the fluid is static and only the radiation 
solver was called by RAMSES, bypassing the hydrodynamic Go¬ 
dunov solver. Both setups used a 2048^ mesh. 

The strong scaling results are displayed in Eig. 5 a. We can 
see that as we go beyond the 12-core limit of the Occigen pro¬ 
cessors, the speedups drop below the ideal curve (grey area), as 
communications begin to take longer to complete. Our radiative 
transfer method appears to perform slightly better than the native 
MHD scheme in RAMSES, as the coupling between hydrodynam¬ 
ics and radiation is ignored (the Godunov solver is not used in 
the RHD simulations). In the weak scaling RHD runs, when the 
size of the problem is doubled with the number of CPUs, it takes 
the implicit BiCGSTAB solver more iterations to converge, as 
there is a stronger propagation of round off errors originating 
from the calls to the MPI_ALLREDUCE routine when more CPUs 
are used. To ensure a fair comparison, we forced all simula¬ 
tions to execute the same number of iterations (500) for every 
timestep, chosen as the maximum observed number of iterations 
in the 1024-core run. All simulations also performed the same 
number of timesteps (100) of a fixed length in time (At = 10“^^ 
s), and each CPU held a grid of 128^ cells. The weak MHD sim¬ 
ulations were run for 300 timesteps, with each CPU processing 
a grid of 512^ cells'^. The results in Pig. 5b show that the weak 
scaling performance of the RHD solver is below the native MHD 
solver. The iterative implicit solver suffers from heavy commu¬ 
nication operations to compute residuals and scalar quantities 
which need to be performed at each iteration for each timestep, 
while the MHD solver only requires one communication opera¬ 
tion per timestep. We believe that a weak scaling efficiency of 
60% remains however acceptable for our purposes. 

4.2. Group scaling 

Our final performance test was to assess the scaling of our multi¬ 
group algorithm for a given problem when the number of fre¬ 
quency groups Ng is increased. The results of the simulation time 
divided by the total number of BiCGSTAB iterations for a ID 
multigroup diffusion test performed on a single Intel® Xeon® 
E5620 (2.40 GHz) CPU core on a local HP-Z800 workstation 
are shown in Pig. 5c (circles). The algorithm appears to scale 
with Ng, which is expected from the double sum over Ng in the 
term on the third line of equation (5). We carried out a sec¬ 
ond group scaling study, this time running the sub-critical ra¬ 
diative shock test from section 3.4 (although with a lower res¬ 
olution; only 3 levels of refinement were used^), where the ra¬ 
diative transfer is fully coupled to the hydrodynamics. The re¬ 
sults are displayed in Pig. 5c (squares), and the behaviour is very 
similar to the diffusion-only solver. The radiative transfer step 

^ https: //WWW .cines.fr/calcul/materiels-et-logiciels/ 
occigen/ 

^ These resolutions and number of timesteps were chosen to that both 
RHD and MHD simulations run for approximately the same amount of 
time. 

^ This explains a smaller time per iteration for the radiative shock sim¬ 
ulations than for the diffusion runs. 


Article number, page 6 of 10 





















M. Gonzalez et aL: Multigroup radiation hydrodynamics with adaptive mesh refinement 


completely dominates over the hydrodynamic step^, in terms of 
computational cost, and it is thus not surprising to see the same 
scaling for a RHD run than for a calculation which only calls the 
radiation solver. 


5. Application to star formation 

In this final section, we apply the multigroup formalism to a 
simulation of the collapse of a gravitationally unstable dense 
cloud core, which eventually forms a protostar in its centre. The 
collapsing material is initially optically thin and all the energy 
gained from compressional heating is transported away by the 
escaping radiation, which causes the cloud to collapse isother- 
mally. As the optical depth of the cloud rises, the cooling is no 
longer effective and the system starts heating up, taking the core 
collapse through its adiabatic phase. A hydrostatic body (also 
known as the first Larson core; Larson 1969) approximately 
10 AU in size is formed and continues to accrete material from 
the surrounding envelope and accretion disk; this first core will 
eventually form a young star, after a second phase of collapse 
triggered by the dissociation of H 2 molecules (see Masunaga & 
Inutsuka 2000, for instance). In this preliminary astrophysical 
application, we shall however focus on the properties of the first 
Larson core for the sake of simplicity. 

We adopt initial conditions similar to those in Commer 9 on 
et al. (2010), who follow Boss & Bodenheimer (1979). A mag¬ 
netised uniform-density sphere of molecular gas, rotating about 
the z-axis with solid body rotation, is placed in a surrounding 
medium a hundred times less dense. The gas and radiation tem¬ 
peratures are 10 K everywhere. The prestellar core mass has 
a mass of 1 M©, a radius Rq = 2500 AU, and a ratio of ro¬ 
tational over gravitational energy of 0.03. To favor fragmen¬ 
tation, we use an m = 2 azimuthal density perturbation with 
an amplitude of 10%. The magnetic field is initially parallel 
to the rotation axis. The strength of the magnetic field is ex¬ 
pressed in terms of the mass-to-flux to critical mass-to-flux ratio 
fi = (M/0)/(M/0)c = 5 (Mouschovias & Spitzer 1976). The 
field strength is invariant along the z direction, and it is 100 
times stronger in a cylinder of radius Rq (with the dense core 
in its centre) than in the surrounding medium^. We used a gas 
equation of state modeling a simple mixture of 73% hydrogen 
and 27% helium (in mass), which takes into account the effects 
of rotational (for H 2 which is the dominant form of hydrogen for 
temperatures below 2000 K) and vibrational degrees of freedom. 
The frequency-dependent dust and gas opacities were taken from 
Vaytet et al. (2013a), assuming a 1% dust content. We used the 
ideal MHD solver of RAMSES, and the grid refinement criterion 
was based on the Jeans mass, ensuring the Jeans length was al¬ 
ways sampled by a minimum of 12 cells. The coarse grid had 
a resolution of 32^, and 11 levels of AMR were enabled, result¬ 
ing in a maximum resolution of 0.15 AU at the finest level. The 
Minerbo flux limiter (Minerbo 1978) was used for these simula¬ 
tions. 

We performed two simulations; one under the grey approx¬ 
imation and a second using 20 frequency groups. The first and 
last groups spanned the frequency ranges (in Hz) [0 ^ 5 x 10^^] 
and [1.3 X 10^^ 00], respectively. The remaining 18 groups 


^ The radiation solver takes up to 90% of the computation time during 
one timestep. 

^ This was chosen to try and reproduce the dragging-in of field lines 
that would have happened in the formation of the dense core (see Gillis 
et al. 1974, for example), while also retaining in the simplest manner 
the divergence-free condition for the MHD. 
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Fig. 6. (a) Temperature as a function of density in the grey and multi¬ 
group simulations at a time t = 24,265 yr. The blue colour represents 
regions where the grey simulation either dominates (in terms of mass 
contained within the figure pixels) over the multigroup simulation, or 
where there is no multigroup data. Likewise, the red codes for the re¬ 
gions of the diagram where the multigroup run prevails. The white areas 
are where both simulations yield identical results. The black contour 
line delineates the region where data are present, (b) Same as for (a) but 
showing the strength of the magnetic field as a function of density. 


were evenly (logarithmically) distributed between 5 x 10^^ and 
1.3 X 10^"^ Hz. The results are shown in Figs 6 and 7. The gas 
temperature as a function of density for all the cells in the com¬ 
putational domain is shown in Fig. 6a, where the blue colour rep¬ 
resents regions where the grey simulation dominates (in terms of 
mass contained within the plot pixels) over the multigroup sim¬ 
ulation, and the red shows where multigroup data prevail. While 
the two simulations show similar results overall, there are several 
differences we wish to point out. For relatively low densities in 
the range 5 x 10“^^ < p < 3 x 10“^^ g cm“^, the multigroup 
run is hotter than the grey simulation. This is also visible in 
the temperature maps of Fig. 7, where the gas is hotter in the 
20-group run for r > 20 AU, this being most obvious in panel 
(c). It appears that the radiation transport from the central core 
to the surrounding envelope is more efficient in the multigroup 
case. More energy has left the core, rending it colder than in 
the grey case, while more energy has been deposited in the thus 
warmer outer envelope. Higher temperatures are also observed 
in the bipolar outflow, along the vertical axis, relatively close to 
the first core. Kuiper et al. (2011) found that using a frequency- 
dependent scheme in simulations of high-mass stars could en¬ 
hance radiation pressure in the polar direction compared to the 
grey simulations of Krumholz et al. (2009), producing much sta- 
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Fig. 7. (a) Temperature map (colours) at a time t - 24,282 yr in a grey simulation of the collapse of a 1 M© cloud. The data are presented as a 

function of radius and height over the mid-plane; the data have been averaged around the azimuthal direction. The black lines are logarithmically 
spaced density contours, in the range 10"^^ < P < 10”^ g cm“^ (two contours per order of magnitude). The arrows represent the velocity vector 
field in the r-z plane, (b) Same as for (a) but using 20 frequency groups, (c) Map of the logarithm of the ratio of the multigroup temperature 
over the grey temperature. The red colour indicates where the gas temperature in the multigroup simulation is higher than for the grey run, and 
vice-versa for the blue colour. 


bier outflows. This could be similar to what we are observing 
here, although this requires further study. 

Conversely, the strength of the magnetic fleld does not 
change signiflcantly when using a multigroup model. In fact, for 
densities below 10“^^ g cm“^, the two runs are virtually iden¬ 
tical (see Fig. 6b). In the ideal MHD limit, the magnetic fleld 
is not directly related to the thermal properties of the gas, and 
it is therefore not surprising that the multigroup formalism has 
little impact. However, we plan in the future to study star for¬ 
mation using a non-ideal description of MHD (Masson et al. 
2012), where magnetic and thermal interactions are twofold. 
First, magnetic diffusion contributes additional gas heating from 
ion-neutral frictions and Joule heating. Second, the chemical 
properties of the gas, which strongly depend on temperature, 
also impact the magnetic resistivities, which govern the diffu¬ 
sion processes. We will investigate this in detail in a forthcoming 
study. 

The use of multigroup RHD may not yield signiflcantly dif¬ 
ferent results in the early stages of molecular cloud collapse, but 
it does provide a wealth of physical information in the system. 
Channel maps such as the ones presented in Fig. 8 or spectral en¬ 
ergy distributions (SEDs; Fig. 9) are directly available from the 
simulation data, without requiring any post-processing software. 
The channel maps show a peak intensity around a wavelength of 
100 yum, close to the synthetic observations of Commer 9 on et al. 
(2012). The SEDs interestingly show departures from a black 
body spectrum, this being most obvious at the low-frequency 
end, close to the protostar (20 AU; Fig. 9b). We do not wish here 


to carry out a detailed study on the effects of multi-frequency ra¬ 
diative transfer on the structures of protostars, we simply wish to 
illustrate the power of the method. We leave the detailed work on 
collapsing objects for a future paper, as this is first and foremost 
a methodology focused article. 


6. Conclusions and future work 

We have implemented a method for multigroup flux-limited dif¬ 
fusion in the RAMSES AMR code for astrophysical fluid dynam¬ 
ics. The method is based on the time-implicit grey ELD solver of 
Commer 9 on et al. (2011b), and uses the adaptive time-stepping 
(in which each level is able to evolve with its own timestep 
using a subcycling procedure) strategy of Commer 9 on et al. 
(2014). The multigroup method allows the discretisation of the 
frequency domain to any desired resolution, enabling us to take 
into account the frequency dependence of emission and absorp¬ 
tion coefficients. The radiative energy density in the frequency 
groups are all coupled together through the matter temperature 
and terms reproducing Doppler shift effects when velocity gra¬ 
dients are present in the fluid. A consequence of this coupling is 
the apparition of non-symmetric terms in the matrix we have to 
invert in our implicit time-stepping procedure. We therefore had 
to abandon the original conjugate gradient algorithm of Com- 
mer 9 on et al. (201 lb) for a bi-conjugate gradient iterative solver. 
A more evolved BiCGSTAB solver was preferred for its greater 
stability compared to a raw bi-conjugate algorithm. 
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Fig. 9. Spectral energy distributions extracted in a cell located 
2000 AU (a) and 20 AU (b) from the protostar. The black solid lines 
represent the energy inside the 20 frequency groups. They are com¬ 
pared to a black body distribution (dashed red) which would have the 
same total energy. 


Fig. 8. Radiative energy density map in each frequency group. The 
white contours have 10 levels logarithmically spaced between the mini¬ 
mum and maximum values of each map. The frequencies of each group 
are indicated at the bottom of each map. The maps represent the same 
region as in Fig. 7. 


The method was fully tested against standard radiation dif¬ 
fusion, frequency-dependent, and full radiation hydrodynamics 
tests. It performed extremely well in all of these tests, and its 
scaling performance was also found to be very satisfactory. 

The multigroup formalism was finally applied to a simula¬ 
tion of the gravitational collapse of a dense molecular cloud core 
in the context of star formation. The method has revealed differ¬ 
ences between grey and frequency-dependent simulations, but 
more importantly uncovered departures from a black-body radi¬ 
ation distribution. We also illustrated the wealth of information 
the method brings to astrophysical studies, with the ability to 
directly produce channel maps and SEDs. We will carry out a 
much more thorough study of the effects of multigroup radiative 
transfer on the structures of protostars and proto-planetary discs, 
as well as their observable quantities, as part of a much wider 
parameter space study in a forthcoming paper. 
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